Transient and Equilibrium Synchronization in Complex Neuronal Networks 
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Transient and equilibrium synchronizations in complex neuronal networks as a consequence of dy- 
namics induced by having sources placed at specific neurons are investigated. The basic integrate- 
and-fire neuron is adopted, and the dynamics is estimated computationally so as to obtain the 
activation at each node along each instant of time. In the transient case, the dynamics is imple- 
mented so as to conserve the total activation entering the system. In our equilibrium investigations, 
the internally stored activation is limited to the value of the respective threshold. The synchro- 
nization of the activation of the network is then quantified in terms of its normalized entropy. The 
equilibrium investigations involve the application of a number of complementary characterization 
methods, including spectra and Principal Component Analysis, as well as of an equivalent model 
capable of reproducing both the transient and equilibrium dynamics. The potential of such concepts 
and measurements is explored with respect to several theoretical models, as well as for the neuronal 
network of C. elegans. A series of interesting results are obtained and discussed, including the fact 
that all models led to a transient period of synchronization, whose specific features depend on the 
topological structures of the networks. The investigations of the equilibrium dynamics revealed a 
series of remarkable insights, including the relationship between spiking oscillations and the hier- 
archical structure of the networks and the identification of twin correlation patterns between node 
degree and total activation, implying that hubs of connectivity are also hubs of integrate-and-fire 
activation. 

PACS numbers: 05.10.-a, 05.40Fb, 89.70.Hj, 89.75.k 



'Measure what is measurable, 
what is not so. ' (G. Galilei) 



and make measurable 



I. INTRODUCTION 

Few research areas have received so much attention 
as neuronal networks and complex networks. Investiga- 
tions at the intersection of these two areas are particu- 
larly promising because they naturally integrate the em- 
phasis on structure typical of complex networks and the 
dynamics systems features of neuronal networks. Such a 
multidisciplinary area is henceforth called complex neu- 
ronal networks. Despite the many possibilities for cross- 
fertilization between these two areas, relatively few re- 
lated works have been reported. Probably the first in- 
vestigations in neuronal complex networks were reported 
in ■ Other related works include Another is- 

sue of special interest regarding the dynamics of systems 
underlain by complex connectivity regards their synchro- 
nization. Synchronization plays an important role in sev- 
eral natural situations, including heartbeats and overall 
brain activity. Several works addressing synchronization 
in complex networks have been reported in the litera- 

ture (e.^[lj7|, I, ©, M M E El M EE EE Ei 

[Til . [T9L 1^ K2lf). The investigation of synchronization 
in neuroscience and neuronal networks has a long his- 
tory (e.g. [H, [H, H3, More recently, growing at- 
tention has been placed on synchronization of complex 
neuronal networks (e.g. [26|, [27], HE However, a 
great deal of such investigations has had their attention 
concentrated on models of individual neuronal activation 



(e.g. Hodgkin-Huxley or Kuramoto) which are more so- 
phisticated that the integrate-and-fire individual neuron 
model. Though many interesting results have been pro- 
duced by such works, it is felt that a more systematic 
exploration of the relationship between synchronization 
and the different types of topologies in complex networks 
can be achieved by using simpler models for the neuronal 
dynamics, such as the integrate-and-fire approach, which 
allows faster computational simulations 

The current article addresses the study of synchroniza- 
tion during the transient activation of complex networks 
with different topological features involving neurons of 
the integrate-and-fire type under conservation of the in- 
coming activity. Although intrinsically related to neu- 
ronal networks, such results extend immediately to other 
types of systems involving integration-and-fire dynam- 
ics, especially production systems in which the output 
of a node requires the integration of specific inputs in 
order to yield a product (e.g. production of molecules, 
industrial production, or computational calculations in 
parallel systems). In this respect, the delay implemented 
by the combination of the integrator and non-linear ele- 
ment would act as a mechanism for favoring the chances 
of receiving all required inputs before the node can act 
on them. 

Specific nodes are chosen as source of activation, and 
their effect in activating other neurons is sound and ob- 
jectively quantified in terms of the instantaneous syn- 
chronization, a measurement defined in terms of the nor- 
malized entropy of the node activations at a given in- 
stant of time. The normalized instantaneous synchro- 
nization (NIS) is related to the entropy of the system 
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activation (e.g. @, H3, HH), but also takes into account 
the number of nodes already activated. Such an analysis 
of the overall synchronization of the network is performed 
with respect to 6 theoretical models of networks (Erdos- 
Renyi, Barabasi- Albert, Watts-Strogatz, a geographical 
network, as well as two knitted networks). The neuronal 
network of C. elegans [f| is also investigated. The ob- 
tained results indicate that the overall synchronization 
tends to present a peak along the medium term, collaps- 
ing subsequently. In addition, each type of network im- 
plied specific features regarding the network activation 
and instantaneous synchronization. 

The patterns of regular spiking arising at the equi- 
librium as a consequence of synchronization of the indi- 
vidual firings are also investigated by using a series of 
characterization methodologies as well as an equivalent 
model. 

This article starts by summarizing the main concepts 
in complex networks and neuronal networks and pro- 
ceeds by describing the normalized instantaneous syn- 
chronization concept, which is subsequently applied for 
the characterization of the synchronization in 6 networks 
of distinct kinds. The equilibrium investigations are 
reported subsequently with respect to Erdo-Renyi and 
Watts-Strogatz Networks. 



II. BASIC CONCEPTS 

This section summarizes the main concepts in network 
representation, measurement, random walks, as well as 
the six network models assumed in the present article. 



A. Complex Networks Basics and Models 

A weighted, directed network T can be fully repre- 
sented in terms of its weight matrix W . Each edge ex- 
tending from node i to node j, with associated weight v, 
implies W(j, i) = v. The absence of connection between 
nodes i and j implies W(j, i) — 0. The out-degree of a 
node i, henceforth expressed as k out , corresponds to the 
number of outgoing edges of that node. The out-strength 
of a node i, s ou t, is given by the sum of the respective 
weights of all outgoing edges. Similar definitions hold for 
the in-degree and in-strength. 

Six models of complex networks are considered 
in the present article: Erdos-Renyi (ER), Barabasi- 
Albert (BA), Watts-Strogatz (WS), a geographical model 
(GG) [HI, [H, HH, as well as the path-regular network 
(PN) and path-transformed model (PA) [H HEJp, [H[ . 
The ER network (see also related works by [39|]) was 
grown by taking each possible edge with constant prob- 
ability, the BA structure was grown by using the tra- 
ditional preferential attachment scheme [32| . and the 
WS was derived from a linear regular network with 
rewiring rate 0.1 0,Q. The geographical network (GG) 



is obtained by distributing the nodes through a two- 
dimensional space and connecting all pairs of nodes which 
have distance smaller than a fixed threshold. The two 
knitted networks are the path-regular (PN) and path- 
transformed B A (PA) networks [36|, l37| . Both these net- 
works are formed by paths. The PN network is grown 
by incorporating paths involving all network nodes, be- 
ing inten sely reg ular regarding several of its topological 
features [33, [38|]. The PA network can be obtained by 
transforming (from stars to paths) a BA network with the 
same number of nodes. All networks in this work have 
similar number of nodes and average degree. Only the 
largest connected component has been taken into account 
for each network. However, because of the relatively high 
average degree adotped in this work ((k) = 6), most of 
the nodes end up belonging to the largest component. 



B. Hierarchical Organization of Complex Networks 

Given a complex network and a reference node i, its 
hierarchical organization can be obtained by flooding the 
network from the reference node [40| . The hierarchical 
organization includes the concentric levels (or hierar- 
chical levels) of the network, namely the levels contain- 
ing the nodes which are at successive shortest path dis- 
tances from the reference node i. So, the first concentric 
level incorporates the original nodes which are at short- 
est path distance 1 from i (i.e. they are the immediate 
neighbors of i). The second concentric level includes the 
nodes which are at shortest path distance 2 from i, and so 
on. The connections between the reference node and the 
nodes at the h— level can be understood as implementing 
virtual links [4l| . Once the hierarchical organization of 
a network has been obtained with respect to a specific 
node chosen as the reference, a series of measurements 
can be calculated S3, |43, including the hier- 
archical number of nodes n^(i), the hierarchical degree 
kh{i) and the intra-ring degree a^(i). The hierarchi- 
cal number of nodes corresponds to the number of nodes 
within each respective concentric level. The hierarchical 
degree is equal to the number of edges from level h to 
level h + 1. The intra-ring (or intra-level) degree corre- 
sponds to the number of edges established within level 
h. 

With the extension of the hierarchical organization to 
complex networks with asymmetric connections [45) . the 
hierarchical degree needs to be split into hierarchical 
indegree and hierarchical outdegree. Therefore, the hier- 
archical indegree of level h with respect to node i, hence- 
forth abbreviated as kih(i), is equal to the number of 
edges received by level h from level h—1. The hierarchi- 
cal outdegree of level h is the number of edges sent from 
that level to level h + 1. 
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C. Instantaneous synchronization 

Let £ be a dynamic system implemented over a com- 
plex network involving TV nodes. The activation of each 
node i at each time t is henceforth represented as A(i, t). 
For simplicity, and without loss of generality, such ac- 
tivations can be normalized so that they become a sta- 
tistical distribution. This can be done by defining the 
probability of activation of node i at time instant t to be 
a(i,t) = A(i, t)/y^ ^_, A(i, t). The entropy of all such ac- 
tivations at t (e.g. [3Q.l3ll|) can not be immediately given 
as 



e(f) = -^A(i,t)log(A(i,t)) 



(1) 



Note that the maximum value of the entropy, corre- 
sponding to log(N), is achieved when all nodes have the 
same probability of activation 1/N, i.e. the total acti- 
vation is the most uniformly distributed amongst all the 
nodes. 

Because the activation probability can be understood 
as a mean frequency of activation (e.g. (20|). it is inter- 
esting to consider the instantaneous synchronization of 
the system. A possibility is to use the following expres- 
sion 1201 



Na 

N 



a(t) 



N a (t) log(N) - e(t) 



N 



log(N) 



(3) 



where N a (t) is the number of nodes with non-zero ac- 
tivity at time t. Now, we will only have maximum syn- 
chronization £(t) = 1 when N a (t) = N and all nodes 
have the same activation probability. This measure- 
ment, henceforth called normalized instantaneous syn- 
chronization (NIS), is adopted throughout this work. Go- 
ing back to the previous example, we now have (,(t) ss 
(0.85) (0.2) « 0.017, which provides a more reasonable 
quantification of the instantaneous synchronization con- 
sidering the whole network. 



D. Spectral Characterization 

Given a time series or signal s(t), it is often quite dif- 
ficult to identify its periodical components. By reinforc- 
ing the intrinsic periodicities along the signal, its auto- 
correlation allows a more effective means for inferring 
its constituent oscillations. In this work we resource to 
the power spectrum P(f) of the signal s(t), which cor- 
responds to the squared magnitude of the Fourier trans- 
form of the autocorrelation of s(t) [46[. More specifically, 
given the signal s(t), its Fourier transform is defined as 



a{t) 



log(N) - e(f) 
log(N) 



(2) 



Note that < a(t) < 1, with the maximum en- 
tropy leading to null synchronization and minimal en- 
tropy leading to maximum synchronization. The nodes 
with null activation are not considered in the calculation 
of the entropy used for the synchronization because they 
are not really participating to the overall activation dy- 
namics. 

Though such a definition properly reflects the rela- 
tionship between activation entropy and instantaneous 
synchronization, it does not take into account the fact 
that, especially during the transient period of time (but 
sometimes also at steady state), some nodes will not be 
active. Let us illustrate this problem through the fol- 
lowing example. Let a dynamical system with N = 100 
nodes have only 2 nodes activated at time t, e.g. be- 
cause we are in the transient period and activation has 
not yet reached the other nodes. The probability acti- 
vation therefore will be a(t) = 0.5, yielding e(t) fx 0.69 
and a = (log(lOO) - e(t))/ log (100) w 0.85, indicating a 
high level of synchronization between the two nodes. Al- 
though this is really the case when only the two nodes are 
considered, the system actually involves other 98 nodes 
which are at zero activation. In order to better express 
the overall synchronization considering all the N nodes, 
we adopt henceforth the following alternative definition 
of the instantaneous synchronization of the dynamics: 



S(f) = / s(t)exp(-i2irft)dt 

J t— — oo 

The autocorrelation of s(i) is 



a(r) = / s(t)s(r - t)dt 

J t— — oo 

Therefore, the Fourier transform of the autocorrelation 
function is immediately given by the correlation theorem 
as 



A(f) = S(f)S(fY 



where l S(f)*' is the conjugate of S. The magnitude of 
the spectrum of s(t) is 



l^(/)l = Vs(f)s(f)*dt 



The power spectrum of s(t) is defined as being equal 
to the squared magnitude of the spectrum of s(t), i.e. 



p(f) = s(f)s(fy 
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The power spectrum allows the identification of auto- 
correlations (especially oscillations) of the original signal. 
More specifically, presence of peaks at specific frequen- 
cies / in the power spectrum indicate the presence of 
respective oscillations with frequency / in the original 
signal s(t). All power spectra in this work are estimated 
by using the Fast Fourier Transform (FFT). 



E. Multivariate Statistical Methods 

The spikes produced along time for each neuron in the 
complex neuronal networks can be understood as pat- 
terns, which can be compared and classified by using mul- 
tivariate statistics [46|, |47], 5H and/or pattern recognition 
methods [46l. |48| . In the present work, we apply the Prin- 
cipal Component Analysis (PCA) methodology [o5l litj ] 
in order to decorrelate the spike patterns [49( and to 
obtain more significant clusters in respective projections 
of the original patterns. 

Given a set of M spike patterns along H time steps, a 
total of H measurements corresponding to the presence 
of a spike at each time can be obtained. Let us define 
the matrix U so that each of its rows corresponds to a 
train of spikes, so that U has dimension M x H. Let ft 
be the 1 x H vector containing the average number of 
spikes at each time step h = 1, 2, . . . , H. Now, define the 
new matrix F as 

F = U - ones(M,l)fl 

where onesiM, 1) is a M x 1 vector of ones. The covari- 
ance matrix, taking into account all pairwise covariances, 
is immediately given as 



C = 



1 



-FF 



M -I 

Observe that C is symmetric. Let Aj be the eigen- 
values of C sorted in descending order, with respective 
eigenvectors wj. The Karhunen-Loeve transform of the 
original measurements is the stochastic linear transfor- 
mation implemented by the matrix G given as follows 



G 



l>2 



(4) 



The PCA method involves transforming the original 
measurements as V — GU T with m << H, as allowed by 
the high redundancy normally found along and between 
signals. It can be shown that the resulting new measure- 
ments V, which correspond to linear combinations of the 
original measurements, are completely decorrelated by 
the PCA methodology, therefore maximizing the varia- 
tion of the data long the first new variables. 



III. MODELING AND SIMULATION 

In this work, we focus attention on the two following 
specific features: (i) integrate-and-fire dynamics at each 
node; and (ii) conservation of incoming activation. Fig- 
ure [1] shows the basic node adopted henceforth, which 
corresponds to a simple integrate-and-fire neuron. Each 
such node i includes n(i) inputs and m(i) outputs. The 
input activity is integrated until its value reaches the 
threshold T(i) (hard limit non-linearity is adopted in this 
work), in which case the neuron fires. 




FIG. 1: The integrate-and-fire type of neuron adopted in this 
work. 



The integrator and non-linear stages are henceforth un- 
derstood to constitute the soma (or body) of the neuron. 
The above type of neuron can be immediately modified 
in order to allow the conservation of the total activa- 
tion which has already entered the system (the system 
may involve sources of activation). First, once the neu- 
ron fires, all its store activation is liberated through the 
outgoing edge (i.e. axons). Observe that the total liber- 
ated activation will be necessarily equal to the threshold 
T(i), and that the activation will be divided amongst 
the outgoing edges, so that each receive a fraction of ac- 
tivation equal to y^j = l/k ou t(i) (so, k ut{i) — m(i)). 
Second, the total activation received by the neuron is 
stored internally until the neuron fires. After each spike 
and transfer of the currently stored activation, the state 
of the neuron is cleared. These two simple modifications 
allow the activation which has already entered the sys- 
tem to be completely conserved along time. It should be 
observed that the division of activation amongst the out- 
going connections is not biologically reasonable because 
the action potential in neurons involves signal reinforce- 
ment (the axon is active) and is known to produce spikes 
with similar intensities in most synapses. However, a 
biologically realistic network can be easily set up in or- 
der to reproduce the conservative dynamics by making 
the synapses at which the axons y^j terminate to have 
weights l/k out (i). 

In the case of the equilibrium analyses, the activation 
inside each memory S(i) is limited to L(i) at all times. 
This implies that the portion of the activation received by 
a node which exceeds L(i) is discarded, therefore under- 
mining the conservation of the activation received from 
the source node after the neurons start firing. 
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The complex neuronal networks considered in this ar- 
ticle consist of a representative sample of each of the 
ER, BA, WS, GG, PN and PA models. The dynamics 
is implemented by considering each of the nodes in these 
networks to be a conservative integrate-and-fire neuron 
as discussed above. All neurons are assumed to have 
the same threshold T = 1. The activation of the net- 
work is implemented by assigning a source to a specified 
node, which therefore acts as a source of constant acti- 
vation with intensity 1. As time passes, such an activa- 
tion is distributed to the other nodes in the networks. 
The total of activated neurons at any time t is N a (t). 
In order to ensure activation conservation, the weight 
of each connection from node i to node j is defined as 
w(j,i) = l/k ou t(i)- For the sake of simplicity, each of 
the undirected edges yielded by the 6 considered net- 
work models are dissociated into one dendrite and one 
axon, so that the out-degree becomes identical to the in- 
degree. Less-symmetric configurations can be considered 
futurely. The C. elegans network is kept directed in our 
simulations. 

The dynamics of such complex neuronal networks has 
been investigated with respect to: (a) the evolution of 
the activity of all nodes along time, represented in dia- 
grams which are henceforth referred to as activograms; 
(b) the distribution of activated and non-activated nodes 
along time; (c) the distribution of the spikes produced 
by all neurons along time (spikegram); and (d) the evo- 
lution of the normalized instantaneous synchronization 
along time. The maximum NIS values obtained while 
considering the source at every node, as well as the time 
at which such values occur, are also considered in this 
work. 



IV. TRANSIENT SYNCHRONIZATION: 
RESULTS AND DISCUSSION 

Simulations on the theoretical models were performed 
considering n — 100 and (k) = 6. The largest connected 
component in the neuronal network of C. elegans con- 
tained 239 neurons. For both the theoretical and real- 
world networks, each of the nodes was considered as a 
source of activation with intensity 1. 

Figures[2]to[7]show the patterns of activation of all neu- 
rons for each time t = 1,2,..., 100 obtained for each of 
the 6 networks by having the source of activation placed 
at node 50 (with a few exceptions, similar patterns were 
identified for the source at other nodes). More specif- 
ically, each of these figures show the activogram, i.e. 
the activation at each neuron along time (a), the active 
(white) and non-activate (black) neurons along time; and 
(c) the spikes produced by each neuron along time. The 
activations, instead of the normalized probabilities of ac- 
tivation, are shown in the activograms in all figures for 
the sake of better visualization. 

A series of interesting results and interpretations can 
be identified from these figures. As expected, all acti- 



vations tended to spread progressively from the source 
node 50 as time passes, with the rate of spikes increasing 
steadily with time. However, quite distinct patterns of 
activation have been observed for each of the considered 
networks. In the case of ER (Fig. [2]), for instance, a rea- 
sonably uniform distribution of activation along time was 
obtained, with most nodes engaging into activity for the 
first time at similar instants (after approximately 10 or 20 
steps, see Figure [2jb)). As shown in Fig.[3l quite a differ- 
ent dynamics of activation has been obtained for the BA 
network. Because of the presence of hubs, several nodes 
are activated relatively soon (less than 10 steps), while 
some nodes are only recruited much later. In addition, 
the activity along time tends to concentrate in the hubs 
at the left-hand side of the image. It would be particu- 
larly interesting to verify whether the rate of individual 
activations follow a power law. Yet another pattern of 
activations has been obtained for the WS model (Fig. 0|, 
which has been made clearer by the fact that the orig- 
inal neurons in the one-dimensional regular lattice used 
to derive this network had been sequentialy numbered. 
Now, the activation proceeds gradually through the suc- 
cessive neighbors. Also, once activated, the neurons seem 
to engage in more regular patterns of spiking than those 
obtained for ER and BA. The onset of activation in the 
GG network is peculiar, especially regarding the fact that 
neurons tend to start activity at the most diverse times. 
This is explained by the fact that the GG is not small- 
world, implying the activation to progress along the adja- 
cencies along the network. Therefore, neurons which are 
connected to the source through longer shortest paths 
will engage into activity later. The activation diagrams 
obtained for the PN and PA networks are remarkably 
similar, being characterized by relatively uniform times 
for activation onset. This result is particularly surprising 
because these two networks are know to have markedly 
distinct structures [37], [H[ • 

The normalized entropies of activations for each time 
considering all neurons, with source of activity at neuron 
50, are shown in Figure [5] With a few exceptions, sim- 
ilar patterns were obtained when the activation source 
was placed at other neurons. These curves tend to be 
similar, involving an initial stage with plateaux of rel- 
atively low entropy, followed by more gradual progres- 
sion to higher entropy approaching the maximum limit 
of log(N) fa 4.61 in all cases. Such initial plateaux 
are mainly a consequence of the distribution of activi- 
ties among the axons of the neurons more immediately 
connected to the source, especially in the cases where 
the node associated to the source, or its more immedi- 
ate neighbors, had large out-degree (recall that this im- 
plies the activity to be distributed amongst the outgoing 
edges). Figure[5]shows the normalized instantaneous syn- 
chronizations (NIS) obtained along time considering all 
neurons, with source of activation at node 50. Interest- 
ingly, after exhibiting some plateaux of synchronization 
at the initial time steps, these curves tend to evolve to a 
peak of synchronization and then decrease to near zero 



FIG. 2: Activation diagrams for ER network: (a) probability 
activation of each node along time (a); the active and non- 
activate neurons along time (b); and the spikes produced by 
each neuron along time (c). The source was placed at node 
50. 



FIG. 4: Activation diagrams for WS network: (a) probability 
activation of each node along time (a); the active and non- 
activate neurons along time (b); and the spikes produced by 
each neuron along time (c). The source was placed at node 
50. 






FIG. 3: Activation diagrams for BA network: (a) probabil- 
ity activation of each node along time (a); the active and 
non-activate neurons along time, shown in black (b); and the 
spikes produced by each neuron along time (c). The source 
was placed at node 50. 




FIG. 5: Activation diagrams for GG network: (a) probability 
activation of each node along time (a); the active and non- 
activate neurons along time (b); and the spikes produced by 
each neuron along time (c). The source was placed at node 
50. 



activation (typically after 60 steps). Observe the higher 
values of NIS obtained for the BA network (Fig. [9b) 
and the more gradual decrease of NIS in the GG case 

(Fig. EH). 

A more complete picture of the instantaneous synchro- 
nization of the networks can be obtained by inspecting 
Figure fTOl where the x— axes correspond to the maximum 
NIS obtained along all times, all nodes, and by consid- 
ering the source in any node; while the y— axes show 
the time w at which the maximum NIS was obtained. 



Remarkably, each of the networks led to relatively ho- 
mogeneous maximum synchronization and time at which 
it occurred, defining reasonably dense clusters of points 
in these scatterplots. The times w at which the max- 
imum synchronizations were observed are similar in all 
cases, except the GG network, and tend to be comprised 
within the interval from 1 to 30 steps. In the case of the 
GG, these times extend to 100 steps, reflecting the fact 
that the GG structure does not present the small-world. 
The highest synchronizations were observed for the BA, 
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FIG. 6: Activation diagrams for PN network: (a) probability 
activation of each node along time (a); the active and non- 
activate neurons along time (b); and the spikes produced by 
each neuron along time (c). The source was placed at node 
50. 




FIG. 7: Activation diagrams for PA network: (a) probability 
activation of each node along time (a); the active and non- 
activate neurons along time (b); and the spikes produced by 
each neuron along time (c). The source was placed at node 
50. 



and the lowest for the WS and GG cases. The curves ob- 
tained for the PN and PA cases are similar one another 
and resemble the ER counterpart. 

Figure [Tl] shows the normalized activations a(i,t) at 
each of the nodes i of the considered GG network at time 
100 (A) and 1000 (b), with the source of activation placed 
at node 50. Observe that the pattern of activation around 
the source tended to be similar in these two cases, sug- 
gesting that the activation tends to a steady-state which 
unfolds more strongly around the source. However, most 



of the other activations changed almost completely along 
the time period between 100 and 1000 steps. 

Figure [T2l shows the activation diagrams (a), activated 
nodes (b) and spikes (c) obtained for all the 239 neurons 
of the C. elegans network along the 100 initial time steps, 
with the activation source at node 50. As can be clearly 
observed from Figure [12l most of the nodes are engaged 
in activity after only about 20 time steps. This seems to 
be related to the relatively high average node out-degree 
(14.36) of this network. Interestingly, the activation of 
the neurons seem to undergo an abrupt increase after 
approximately 100 steps, with most neurons presenting 
similar frequencies of activation thereafter. In addition, 
several neurons tended to exhibit similar frequency of 
spikes after such a transition. The peak of activation 
at the latest stages occurs at the maximum hub of this 
network (out-degree 73). 

The evolutions of entropy and NIS along time with ac- 
tivations placed at nodes 1 to 50, shown in Figure IT51 
are particularly diversified, with initial plateaux present- 
ing quite different lengths and heights. Interestingly, the 
lengthier plateaux also tended to yield the highest nor- 
malized instantaneous synchronization. Also noticeable 
are the shorter and lower plateaux which appear from 
time stapes 1 to 20. However, after nearly 60 steps the 
synchronizations collapse. 

The maximum NISs and respective times at which they 
occurred for the activation placed at each of the nodes 
in the C. elegans network are shown in Figure 1141 It can 
be inferred from this figure that the C. elegans presents 
moderate values of maximum normalized instantaneous 
synchronizations, which occurs at an intensity compara- 
ble to that obtained for the ER, PN and PA networks. 
Interestingly, in some cases the maximum synchroniza- 
tion took place as late as at the 67th and 73th steps. 



V. OSCILLATIONS IN UNIFORMLY-RANDOM 
NETWORKS 

In this section we consider a series of complementary 
approaches in order to characterize the transient and 
equilibrium dynamics in two types of integrate-and-firc 
complex neuronal networks — namely ER and WS struc- 
ture. Representative samples of each of theses types of 
networks have been selected with sizes N = 25, 50 and 
100 and average degrees (k) = 10,20,30,40 and 50. 

Figure [TH] shows the spikegrams obtained for the sev- 
eral configurations of ER networks. Observe the sell- 
defined avalanche transitions at the beginning of the 
spikes obtained for all network configurations. In agree- 
ment with previous results [Hoj , the avalanches initiation 
times were found to be quite similar for all network sizes, 
irrespectively of the average degree. The early spikes ob- 
tained for N = 50 and (k) = 10; N = 100 and (k) = 10; 
N = 100 and (k) = 20 and N = 100; (k) = 30 and 
(k) = 50 were a consequence of the existence of commu- 
nities which arise in the ER networks arising from ran- 
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dom fluctuations. Interestingly, the oscillations tended to 
become more regular (more similar frequencies) and syn- 
chronized with the increase of the average degree. The 
routes to regularity and synchronization this parameter 
is increased exhibited groups of neurons producing simi- 
lar spiking patterns. 

The total number of spikes, shown in Figure [TH] re- 
flects the ensemble behavior of the respective complex 
neuronal networks, arising as a consequence of the lin- 
ear superposition of the spikes being produced by each 
neuron (analogous to EEG potentials). The avalanche 
transitions can be clearly identified, taking place after 
nearly 40 steps for N — 25, 40 steps for N — 50 and 
80 steps for N = 100, which confirms the previous study 
reported in [5(j. In several cases, a clear intense peak 
is observed at the avalanche time, which is followed by 
oscillations whose regularity tend to increase with the 
average degree. 

In addition to the total number of spikes, it is also im- 
portant to take into account the respective power spectra, 
which are shown in Figure 1171 Only the second half of 
each signal (number os spikes) has been considered for 
the calculation of the power spectra in order to represent 
the equilibrium regime. So, a total H — 500 time steps 
have been considered for the estimation of each spectrum. 

Relatively rich spectral composition can be observed 
for most cases, suggesting particularly complex behavior 
of the oscillatory components. The progressive elimina- 
tion of frequencies with the increase of the average de- 
gree is evident from Figure [T71 especially for N = 50 



and N = 100. In addition, the peaks of the spectra for 
N = 50 and N = 100 tend to appear near 20 and 10 fre- 
quency units, respectively. At the same time, it should 
be observed that the peaks of the spectra in Figure fTTl 
tend to shift to the right-hand side for larger values of 
the average degree, signaling the increase of the main 
frequencies. The more regular oscillations obtained for 
larger average degrees are related to the fact that the 
higher this parameter, the more regular the degrees of 
the networks become. At the limit, for very large av- 
erage degree, all nodes become connected, implying a 
highly regular structure as far as all possible topological 
features are concerned. 

Unlike in the previous works [H, US HH, HI, Hl| i the 
activation inside each neuron is limited to the respec- 
tive maximum value L(i). In addition to being more 
biologically-realistic, such a choice also allowed particu- 
larly interesting oscillatory dynamics. At the same time, 
the limitation of the internal activation implies that the 
overall activation, constantly received from the source 
node, is no longer guaranteed to be conserved. It is there- 
fore interesting to consider the total activation inside the 
whole complex neuronal networks along time. Such a 
measurement is depicted in Figure [T51 Interestingly, two 
clearly distinct regimes are immediately identified: one 
transient period in which the internal activation increases 
linearly, followed by the steady-state regime character- 
ized by nearly constant activation (even though activa- 
tion is continuously pumped into the system through the 
source node). The transition between these two regimes 
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FIG. 9: The normalized instantaneous synchronization (NIS) along time for each of the considered networks. The source was 
placed at node 50. 



is in most cases characterized by one abrupt surge of ac- 
tivation, related to the main avalanche. Such a result 
suggests that the main avalanche represents watershed 
signaling the transition of the dynamics in the system 
from conservative to dissipative. 

Because the train of spikes produced by each individ- 
ual neuron, shown in Figure [T51 seems to be organized in 
clusters during the route towards more regular and syn- 
chronous firing (i.e. similar trains of spikes are obtained 
for groups of neurons) , it is interesting to investigate such 
a possibility further. We do this with the help of the PC A 
method, which produced the results shown in Figure [T9l 

It is clear from this figure that the increase of the aver- 
age node degree tends to produce denser groups of spik- 
ing patterns. At the same time, interesting clustering 
structures are obtained for small average node degrees, 
incorporating a dense central cluster (e.g. the PCA scat- 
terplots for N = 25 and (k) = 10, N = 50 and (Jfe) = 10, 
and N = 100 and (k) = 10) surrounded by outliers. The 
dense cluster of spiking patterns tend to survive even for 
relatively large average degrees in the case of N = 100. 

The correlations between the total number of spikes 
produced by each neuron during the simulations and 
the respective degrees have also been considered, yield- 
ing some of the most interesting results reported in the 
present work. These scatterplots are shown in Figure |2T)1 
Striking twin-linear correlation patterns have been ob- 
tained throughout. Such scatterplots indicate that a 
node with a given degree can spike at two distinct fre- 
quencies of spikes, with the number of spikes being lin- 



early related to the node degrees. At the same time, 
given one of such rules, the number of spikes (i.e fre- 
quency) tends to be correlated with the degree of the 
original nodes. Such twin correlations are likely to be 
related to clusters of original nodes appearing at differ- 
ent concentric hierarchical levels, as discussed further in 
Section US 



VI. OSCILLATIONS IN SMALL- WORLD 
NETWORKS 

Having obtained a comprehensive characterization of 
the oscillations induced along the equilibrium regime in 
ER networks, it is interesting to shift our attention to 
the steady-state dynamics of WS integrate-and-fire acti- 
vations. Figure |2"T1 shows the spikegrams obtained for the 
several parametric configurations of the WS models. 

Several interesting features are evident from the re- 
sults shown in Figure [2U Foremost are the more grad- 
ual avalanches, which tend to be centralized around the 
source node (node 1) [55| . Longer avalanche initiation 
times are implied by larger network sizes N. Because 
the increase of the average degree implies more degree- 
regular networks, the gradual activation of the neurons 
observed for relatively smaller average node degrees be- 
came sharper for larger average node degrees. As with 
the ER cases, the spiking patterns also tended to become 
more regular and synchronized for larger values of aver- 
age degree. 
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FIG. 10: Scatterplos of the maximum NIS (y— axes) and time at which it manifested itself (x— axes) obtained for each networks 
by considering activations at each possible node. 



Figure [22] depicts the power spectra obtained for each 
of the considered WS configurations. Such spectra, which 
were obtained by considering the equilibrium regime of 
the spikings (i.e. the time steps from 500 to 1000) are re- 
markably similar to those obtained by the ER structures 
(Fig. \T7\ suggesting some possible universality between 
different complex networks types. 

Regarding the total activation inside the system, the 
revealed results (not shown in this work) also included 
the two regimes as in the ER case, but with sublinear 
increase during the transient regime, which was also duly 
followed by the dissipative plateau of overall activation. 

Figure [20] shows the scatterplots of the original node 
degrees and total number of spikes induced by each in- 
dividual neuronal cell. Though these two features also 
tended to be organized in terms of parallel straight 
groups for small values of average degree, the correla- 
tions are now substantially less definite. 



VII. THE EQUIVALENT MODEL FOR 
NON-UNIFORM DEGREES 

Figure [2U illustrates the procedure which had to be 
adopted in order to obtain the a more effective version of 
the equivalent model of complex networks, by taking into 



account non-uniformities of degrees amongst the original 
nodes. First, the hierarchical organization of the origi- 
nal network with reference to a given node i is obtained 
fFig |2~4"]). In this specific case we have 5 concentric lev- 
els (h = 0, 1, . . . , H = 4), each containing 1, 3, 6, 2 and 1 
nodes (i.e. no(i) = 1; ni(3) = 1; 112(1) — 6; 713(1) = 2; 
rii(i) = 1). The nodes within each concentric level which 
have identical degrees are then subsumed into equivalent 
nodes. Thus, nodes 3 and 4 at level h = 1 — all with 
degree equal to 5 — become associated to the equiva- 
lent node B. Nodes 5, 6 and 7, which have degree 4, 
are subsumed by node F. The topology of the resulting 
equivalent model is illustrated in Figure refngiequiv. 

Now, the weights associated to the edges in the equiv- 
alent structure in Figure [M] are determined, similarly as 
done in [451 ] . by means of the following equation, applied 
to each pair of equivalent nodes v and p 

W(v,p) — k(v,p)/d 

where k(v,p) is the number of original edges going from 
the nodes associated to the equivalent node p to the nodes 
associated to the equivalent node v, d = J2 g en(p) Mff>f) 
and f2(p) is the set of equivalent nodes that receive a di- 
rected edge from p. Obseve that W(b, a) expresses the 
weight of the connection from a to b. The weights there- 
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fore obtained for the graph in Figure [24] are given as 
follows 



1/5 1/3 

2/3 1/5 1/2 1/4 

1/3 1/8 

11= 1/5 

2/5 1/4 4/6 

2/3 3/8 1/6 1 

1/2 1/6 



For instance, as the two original nodes associated to 
the equivalent node B, both of which with degree 5, send 
two edges to the equivalent node A, two edges to D, 
four edges to E, while two edges remain inside B, we 
have W(A,B) = 1/5; W(D,B) = 1/5; W{E,B) = 2/5; 
W{B,B) = 1/5. 

The thresholds associated to neuron corresponding to 
each equivalent node are immediately given as the num- 
ber of original nodes associated to that equivalent node. 
Therefore, we have that T(A) = 1; T(B) = 2; T(C) = 1; 
T(D) = 2; T(E) = 1; T{F) = 3; T(G) = 2; and 
T(H) = 1. The limitations to the activations stored into 
each memory are also set as being equal to the respec- 
tive number of nodes within each equivalent node, i.e. 
T(i) — L(i). Observe that the so-obtained equivalent 
network itself correspond to an integrate-and-fire com- 
plex neuronal networks with varying thresholds. 



VIII. EQUIVALENT MODEL PREDICTIONS 

Having characterized the oscillatory integrate-and-firc 
dynamics along the transient and equilibrium regimes in 
terms of several measurements, it is time to consider the 
respective equivalent models. The importance of such an 
approach lies in the fact that these structures involve a 
fraction of the original number of nodes, therefore allow- 
ing the identification and theoretical modeling of more 
definite relationships between structure and dynamics. 
We start by considering the distribution of the equivalent 
nodes and weights, as well as by comparing the predicted 
and real number of spikes and spectra for ER. 

Tables U and |TT] give several features of the topologi- 
cal structure of the equivalent models obtained for each 
of the ER and WS configurations, respectively. These 
features include the total number of original nodes in 
each concentric level (upper line), the hierarchical de- 
grees (second line) , the total number of equivalent nodes 
(in bold, third line), and the number of equivalent nodes 
per concentric level (within brackets, third line). It is 
clear from such measurements that the WS structure im- 
plied a larger number of concentric levels in all cases. On 
the other hand, most of the ER structures implied two or 
three levels (the first of each level corresponding to the 
source node), except for the network with N — 100 and 
(k) = 10, which implied 4 concentric levels. Also, because 



of its greater degree uniformity, the WS model yielded 
substantially fewer equivalent nodes than the ER struc- 
tures. The relatively large number of equivalent nodes 
obtained for each concentric level confirms the fact that 
these two types of supposedly regular networks do actu- 
ally a wide dispersion of node degrees, especially in the 
case of the ER configurations. That is the reason why 
equivalent models which do not take into account such a 
heterogeneity cannot yield accurate predictions. 

Though the partitioning of the concentric levels in 
terms of groups of nodes with identical degrees was crit- 
ical for proper modeling of the oscillatory dynamics, it 
supplied as byproduct a new hierarchical measurement 
corresponding to the number of nodes with specific de- 
grees at each concentric level. Because such features have 
been presently found to be essential for the integrate-and- 
fire dynamics, it is likely that they will also provide valu- 
able resources for the characterization of the topology of 
the networks. 

The weights of the equivalent models obtained for each 
of the ER configurations are shown in Figure [55] The 
equivalent nodes are organized bottom-upwards along 
the columns and rows according to increasing degrees 
of the respectively subsumed original nodes. It is clear 
from these matrices that there are some equivalent nodes 
which receive converging intense weights (the clearer rows 
in the matrices — recall that we adopt the convention 
that the edges extend from the columns towards the 
rows). Such nodes are particularly important for the 
non-linear dynamics because of their greater tendency 
to spike, defining the main avalanches. 

We now present some predictions of the dynamical fea- 
tures as obtained by using the respective equivalent mod- 
els. Figure [26] shows the number of spikes predicted by 
the equivalent model for all the considered ER config- 
urations. By comparing with the respective real num- 
ber of spikes shown in Figure [121 it becomes clear that 
the equivalent model allowed an impressive estimation of 
both the transient and steady-state dynamics in every 
case. 

Similarly impressive predictions of several aspects of 
the dynamics, including spectra, were also obtained for 
the WS configurations but are not shown here. 



IX. EQUILIBRIUM CONFIGURATIONS: 
OVERALL DISCUSSION 



At this point, after having characterized several as- 
pects of the integrate-and-fire dynamics obtained for the 
original networks, as well as the features of the respective 
equivalent models and predictions, it is time to integrate 
all such information in order to obtain a more compre- 
hensive explanation of the origin and properties of the 
oscillations observed for the integrate-and-fire complex 
neuronal networks. 

The two clusters of correlations observed for the ER 
model (Fig. I2T)]) are particularly relevant, suggesting that 
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1, 17, 7 


1, 23, 1 
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24, 
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16: (1) (10) (5) 


5: (1) (3) (1) 


2= (1) (1) 


2= (1) (1) 


2= (1) (1) 
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= 50 


1, 14, 35 


1, 35, 14 


1, 42, 7 


1, 45, 4 


1, 49 






14, 168, 


35, 310, 


42, 246, 


45, 176, 
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22: (1) (13) (8) 


17: (1) (12) (4) 
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2= (1) (1) 
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1, 12, 83, 4 
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36, 859, 


67, 1117, 


68, 1333, 


75, 1358, 






31: (1) (9) (17) (4) 


34: (1) (15) (18) 


33: (1) (19) (13) 


36: (1) (20) (15) 


34: (1) (20) (13) 



TABLE I: Features of the equivalent models obtained for the ER configurations: number of nodes per level (first line); 
hierarchical degrees of each level (second line); and total number of equivalent nodes (bold, third line) and the number of 
equivalent per concentric level (within brackets). 





(k) = 10 


(k) = 20 


(k) = 30 


(k) = 40 


(k) = 50 


N 


= 25 


1, 10, 12, 2 


1, 18, 6 


1, 21, 3 


1, 23, 1 


1, 21, 3 






10, 32, 15, 


18, 84, 


21, 54, 


23, 21, 


21, 57, 






9: (1) (2) (4) (2) 


10: (1) (6) (3) 


9: (1) (5) (3) 


8: (1) (6) (1) 


10: (1) (6) (3) 


N 


= 50 


1, 9, 14, 24, 2 


1, 19, 26, 4 


1, 28, 21 


1, 38, 11 


1, 46, 3 






9, 36, 54, 17, 


19, 131, 62, 


28, 271, 


38, 309, 


46, 114, 






15: (1) (4) (4) (4) (2) 


17: (1) (6) (8) (2) 


16: (1) (8) (7) 


16: (1) (9) (6) 


12: (1) (8) (3) 


N - 


= 100 


1, 11, 34, 38, 16 


1, 21, 57, 21 


1, 27, 55, 1 


1, 42, 57 


1, 50, 49 






11, 57, 118, 91, 


21, 179, 260, 


27, 312, 338, 


42, 577, 


50, 843, 






18: (1) (3) (4) (6) (4) 


21: (1) (7) (9) (4 


23: (1) (8) (9) (5) 


22: (1) (11) (10) 


19: (1) (9) (9) 



TABLE II: Features of the equivalent models obtained for the WS configurations: number of nodes per level (first line); 
hierarchical degrees of each level (second line); and total number of equivalent nodes (bold, third line) and the number of 
equivalent per concentric level (within brackets). 



the spiking frequency tends to increase linearly with the 
indegree of the neurons. Indeed, the higher the indegree, 
the more activation a neuron will receive along time, en- 
hancing its chance of firing. Interestingly, the nodes be- 
longing to each of the two groups straight correlations in 
(Fig. [201) have been found to belong to distinct respective 
hierarchical levels. As the neurons at different concentric 
levels tend to fire with different frequencies, depending 
on their respective number of nodes and connections, the 
two groups of correlations are obtained, with the intra- 
group straight dispersions being accounted by the above 
observed tendency of the spiking frequency to increase 
linearly with the indegrees. 

Despite the clear separation of the hierarchical levels 
provided by the correlation diagrams in Figure 1201 it is 
remarkable that the spiking patterns in Figure [BJ] failed 
completely to clusterize with respect to the hierarchies. 
This interesting fact can be taken as an indication that 
though the nodes belonging to different hierarchical lev- 
els do present a well-defined mean frequency, they are 
highly irregular as far as the interspike times are con- 
cerned. This can indeed be corroborated by the visual 
analysis of the spikegrams in Figure [TS] as well as from 
the rich spectra (in the sense of exhibiting many frequen- 
cies) shown in Figure [T71 With this respect, it would be 



interesting to consider PCA projections taking into ac- 
count smoothed versions of the spiking patterns, which 
would enhance the correlations between those patterns 
and perhaps make more evident the hierarchical clusters. 

Because the considered WS structures implied larger 
number of hierarchical levels (see Tab. [TT|, with respec- 
tive neurons firing at distinct group frequencies, the lin- 
ear relationship between the spiking frequency and the 
indegree became blurred in the respective scatterplots in 
Figure [H 



X. CONCLUDING REMARKS 

A. Transient Synchronization 

The subjects of complex networks, neuronal networks 
and synchronization have special importance in the in- 
vestigation of complex systems and natural phenomena. 
The current work has brought these three issues together 
with respect to transient non-linear dynamics unfolding 
in complex networks with different structures and con- 
servation of activation. More specifically, the normalized 
instantaneous synchronization (NIS) has been proposed 
as a measurement of the instantaneous synchronization of 
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the dynamics among the nodes in the networks at each 
specific time. The activation of the networks was per- 
formed by placing a source of unitary activity at each 
specific node, and the respective dynamics observed and 
characterized in terms of the maximum NIS values, as 
well as the respective times when they occurred. 

The obtained results indicate that the normalized in- 
stantaneous synchronization tends to increase along the 
initial steps and then collapse. Also, the intrinsic topo- 
logical organization of each of the considered types of 
networks was verified to imply markedly diverse patterns 
of activation and maximum NIS. While relatively uniform 
patterns of activation spreading were observed for the ER 
model, groups of hubs tended to concentrate the activity 
in BA networks. Also, the onset of activation of nodes 
was verified to be more uniform for the ER case than for 
BA. The WS structure yielded a pattern of activation 
which tended to spread gradually amongst the neighbors 
of the source. Similar activation and synchronization 
were observed for the two knitted networks, namely PN 
and PA. This is particularly surprising because, though 
both these networks are defined in terms of paths, they 
have completely different degrees of regularity [33, l38j ]. 
The different networks also implied distinct maximum 
NISs, with the BA resulting more synchronized along the 
transient dynamics of activation. Quite diverse times of 
maximum activation were observed for the geographical 
network. Regarding the C. elegans network, it was found 
to exhibit diverse dynamics with respect to the position 
of the source. In addition, most nodes tended to start 
activity with remarkable uniformity during the initial 20 
time steps. The overall activation in this network under- 
went an abrupt increase after nearly 100 time steps. At 
the longest term, the hubs tended to dominate the acti- 
vation dynamics. It is interesting to observe that the ap- 
proaches developed in this work are relevant not only for 
the synchronization studies, but also for the characteri- 
zation of the activation in non-linear systems underlain 
by complex connectivity. 

The perspectives for future investigations are varied. 
Among the possibilities, it would be interesting to con- 
sider other types of activations, e.g. involving sources at 
more than one node or periodic activation instead of the 
constant values used in this work. It would be particu- 
larly interesting to study the combined potential of spe- 
cific sources for defining diverse dynamical features of the 
neuronal activity, especially regarding the facilitation of 
one source of activation by other sources. Another pos- 
sibility is to consider the synchronizability of the rates 
of accesses to a specific node from activity originating 
at several nodes [20(. Such investigation, which is al- 
lowed by the conservation of activity [561 ] , would involve 
the identification of the original source of activations, as 
well as its displacement along the networks, during the 
dissemination of the activation. Other interesting ques- 
tions concern the instant frequency of spiking along time 
for each node, as well as the quantification of correla- 
tions and other types of relationships between the acti- 



vations. Because the activation of most of the considered 
networks tends to undergo an abrupt dissemination af- 
ter an initial transient period, it would be interesting 
to investigate for possible critical dynamics (e.g. phase 
transition). It would also be useful to characterize the 
steady state of activations. Although the concepts and 
methods reported in this article have been considered 
from the specific perspective of neuronal networks, they 
can be immediately extended to investigations of other 
situations such as those involving cortical and biological 
systems, particularly gene activation and protein synthe- 
sis. 



B. Equilibrium Synchronization 

Though relatively simple, the integrate-and-fire model 
yields rich dynamical features, including avalanches, ac- 
tivation confinement inside communities and oscillations. 
Having addressed and explained the first two phenomena 
by using equivalent models [45|, H3| , it is now interesting 
to consider the origins and properties of the oscillations 
observed at both transient and equilibrium regimes in 
integrate-and-fire complex neuronal networks. This con- 
stituted precisely the objective of the present work. Its 
main contributions are reviewed and discussed as follows. 

Characterization of Several Aspects of the Os- 
cillations: The oscillations in the integrate-and-fire 
complex neuronal networks were characterized in terms 
of several measurements and approaches, including the 
visualization of the spikegrams, the total number of 
spikes along time and respective power spectra, total acti- 
vation along time, PCA decorrelation of the spiking pat- 
terns in order to seek for clusterized dynamics features, 
as well as correlations between the total number of spikes 
and degrees of each neuron. Each of these approaches al- 
lowed interesting complementary insights about the os- 
cillatory behavior in uniformly-random (ER) and small- 
world (WS) configurations. Particularly revealing were 
the twin correlations observed for both models, which 
were found to correspond to different frequencies of os- 
cillations taking place at different concentric levels of the 
networks. Such correlations also indicated that the mean 
spiking frequency is linearly related to the indegree of the 
respective neurons, which constutes a particularly rele- 
vant property of the analyzed structures. Analogously 
to the investigation of linear diffusive dynamics reported 
in [54| . this property links structure and oscillatory dy- 
namics, implying that topological hubs will also become 
hubs of activity in the integrate-and-fire dynamics. The 
scatterplots obtained by the PCA projection of the in- 
dividual spiking patterns yielded interesting structures 
which, however, were not in correspondence with the 
twin-correlation partitioniongs. The spectra, rich in fre- 
quencies, confirmed the complex structure of the spik- 
ings. 

Identification of Two Clearly-Defined Regimes: 

The consideration of the total activation along time re- 
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vealed two distinct and well-defined regimes: a conserva- 
tive transient period, followed by the dissipative steady- 
state regime. The critical point separating these two dy- 
namics was found to correspond to the main avalanches. 
The transient evolution was found to be more gradual for 
the WS configurations, which could be indeed expected 
because of the less pronounced (if any) avalanches in this 
type of networks. 

Equivalent Model for Non-Uniform Degrees: 
Following the comprehensive characterization of the os- 
cillatory properties of the integrate-and-fire dynamics 
in complex neuronal networks, an enhanced equivalent 
model was developed which takes into account the diver- 
sity of degrees within each concentric topological level. 
More specifically, the nodes with identical degrees at each 
of the levels are associated to respective equivalent nodes, 
allowing a more detailed model. Such a modification was 
found to be critical for paving the way to impressively ac- 
curate predictions of both the transient and steady-state 
features of the non-linear dynamics. 

Limitation of the Stored Activation: Differently 
from the previous approaches reported in (45l . [Hoi . l5ll . 152 . 
f53j . the activation stored inside the state (memory) of 
each neuron was limited to be at most equal to the respec- 
tive threshold. In addition to being more biologically- 
realistic in which concerns neuronal networks, this choice 
also allowed more interesting oscillatory dynamics. 

Additional Hierarchical Measurements of Net- 
work Topology: In addition to paving the way for pre- 
cise estimations of the dynamics, the enhanced equivalent 
model also yielded as a byproduct new hierarchical mea- 
surements corresponding to the number of nodes with 
specific degrees found at each of the respective concen- 
tric levels. Because such topological features proved to be 
fundamental for modeling the non-linear integrate-and- 
fire dynamics, it is expected that they can also provide 
valuable features for the characterization and classifica- 
tion of complex networks [35| . 

Identification of the Origin and Properties of 
Oscillations: As a consequence of the comprehensive 
characterization and modeling of the integrate-and-fire 
dynamics, a better understanding of the origins and prop- 
erties of the spiking oscillations have been achieved. Of 
special importance are the positive correlations between 
the mean spiking frequency and the neuron indegrees and 
the fact that nodes at distinct concentric levels tend to 
fire with different ensemble frequencies. The combina- 
tion of these two effects yielded the remarkable twin- 
correlation diagrams. In addition, it has been found that 
the oscillations tend to unfold after the main avalanches, 
i.e. along the steady-state of the dynamics, as revealed 
by the total activation in terms of time. 

Several are the possibilities for further related investi- 
gations. In principle, most of the suggestions for future 
work identified in (45|, [5(| HH, [H, [H[ can be immediately 
extended to the present work. More specific possibil- 
ities for futher investigations are identified and briefly 
discussed in the following. 



Further Simplifications by Considering Degree 
Intervals: It would be interesting to consider the sub- 
suming of the nodes in each concentric level not strictly 
under the condition of identical degrees, but by having 
similar degrees (i.e. comprised within specific intervals). 
Such a modification would immediately contribute to 
a further reduction of the overall number of equivalent 
nodes, making the equivalent models more compact. It 
would be interesting to quantify the effect of such a sim- 
plification on the accuracy of the respective predictions 
of the dynamical features. 

Modular Equivalent Model with Non-Uniform 
Degrees: Because of the finer level of representation 
allowed by the currently reported equivalent model, it 
would be potentially useful to extend the previous equiv- 
alent models used to explain and predict avalanches and 
activation confinement inside communities. 

Oscillation Analysis with Action Potentials: 
While the current work, analogously to the previous 
works [111, HH [H, [Hj], considered the activation of 
each spike to be equally distributed amongst the outgo- 
ing edges (i.e. axons), it would be interesting to consider 
the more biologically-realistic hypothesis that the spikes 
have constant amplitude (action potential). Preliminary 
investigations have already indicated that avalanches, ac- 
tivation confinement and oscillations are all present in 
this type of non-linear dynamics. However, it would be 
interesting to perform more systematic related investiga- 
tions. 

Chaos and Chaos Control: The particularly com- 
plex spiking patterns obtained after the main avalanche 
suggest that the dynamics of the considered networks 
may be chaotic. It would be particularly promising to 
apply methods from dynamics systems, such as delay 
diagrams and fractal dimensions, in order to search for 
chaotic behavior in the integrate-and-fire complex neu- 
ronal networks. Another interesting perspective would be 
to apply concepts from chaos control in order to interfere 
with the several remarkable aspects of the investigated 
dynamics. 

Distinct Weights: So far, we have been limited to 
integrate-and-fire complex neuronal networks containing 
identical weights. Indeed, the weights have only been 
used in order to obtain the respective equivalent models. 
In this case, the weights correspond to the proportions of 
the activations (edges) which are sent to each equivalent 
node. It would be interesting to consider complex neu- 
ronal networks incorporating weights, in order to inves- 
tigate the effects of such distributions on the respective 
dynamics. 

Investigation of Oscillations in Other Theoret- 
ical Complex Networks Models: Because of space 
restrictions, the oscillatory behavior of integrate-and-fire 
complex neuronal networks described in this work has 
been restricted to a uniformly-random and a small-world 
theoretical models of complex networks. It would be par- 
ticularly promising to extend such an investigation to 
other important types of network topologies, especially 
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scale free and geographical. The consideration of the 
highly regular knitted structures [37| would also be in- 
teresting. 

Applications to Real- World Problems: Though 
we have so far considered the oscillations in theoretical 
models of complex neuronal networks, the respective con- 
cepts and methods are immediately applicable to real- 
world networks. Of particular interest would be to char- 
acterize and model the neuronal network of C. elegans, 



as well as cortical networks. 
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FIG. 11: The normalized instantaneous synchronizations at 
each node of the GG network considered in this article at 
times 100 (a) and 1000 (b). 
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FIG. 12: The activation diagrams for the 239 nodes in the largest connected component of the C. elegans network. 
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FIG. 13: The entropies and NISs along time for neurons 1 to 
50 in the C. elegans network. 
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FIG. 14: The maximum NIS and respective times obtained 
along the initial 100 steps for the C. elegans network, consid- 
ering the source placed at each of the nodes. 
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FIG. 15: The spikegrams, i.e. the occurrence of spikes for each neuron (identified along the a;— axis) along time (y— axis) 
obtained for the ER configurations. Observe the avalanche transitions during the transient regime, as well as the increase of 
frequency regularity and synchronization with observed for larger values of average degree. 




FIG. 16: The number of spikes along time obtained for the original ER complex network configurations. 
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FIG. 17: The power spectra obtained for the number of spikes produced by the whole original ER networks. 




FIG. 18: The total activation inside the complex neuronal networks is not conserved (recall that activation is being continuously 
fed from the source node) after the main avalanche, but rather reaches a plateau of activation after the transient regime. 
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FIG. 19: The two-dimensional scatterplots obtained by the application of the PCA method over the spike patterns of each 
neuron in the ER configurations. 
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FIG. 20: The relationships between the total number of spikes per time step generated by each neuron and the respective degrees 
for all the considered ER configurations. The number of spikes were counted during the last 500 steps of 1000 steps-long signals. 




FIG. 21: The spikegrams obtained for the WS configurations. 
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FIG. 22: The power spectra obtained for the number of spikes produced by the whole original WS networks. 
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FIG. 23: The relationships between the total number of spikes per time step generated by each neuron and the respective 
degrees, all the considered WS configurations. The number of spikes were counted during 1000 time steps. 
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FIG. 24: The hierarchical organization of a simple network (a) with respect to node 1. A total of 5 concentric levels have 
been obtained. The dotted boxed identify the nodes with identical degrees at each concentric level, which are subsumed by the 
respective equivalent nodes, yielding the equivalent network shown in (b). 
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FIG. 26: The total number of spikes along time t for the ER configurations as predicted by the respective equivalent models. 



